load("Implementation Main Data.Rdata")

#Main Implementation Results - Figure 6 #####
controls<-c("truth_init","prisr_init", "powtran_init",
            "as.factor(incompatibility)",
            "spline","spline2","spline3")

base_form<-as.formula(paste("sig_minor_war_mh~as.factor(bin)*as.factor(amnest_init)+",
                            paste(controls, collapse="+")))
base<-brglm(base_form, data=data)
set.seed(1015)
base_boot<-cluster_bootstrap(base, n=1000, amn="amnest_init", bin="bin")
m1<-int_plot(base, base_boot, amn="amnest_init", bin="bin")
m1[[1]]

pdf("Fig6.pdf", width=7, height=4)
print(m1[[1]])
dev.off()

#Full Main Results Table - Appendix Table 7#####
set.seed(1015)
n<-1000
coef<-matrix(NA, ncol=length(coef(base)))

pop<-unique(data$pam_caseid)
i<-1
model<-NA
while(nrow(na.omit(coef))<n){
  #Make within loop holder for effects
  effects<-matrix(NA, ncol=length(coef(base)))
  sample<-sample(pop, replace = T, size = length(pop))
  map<-sapply(sample, function(x) which(data[,"pam_caseid"]==x))
  boot.sample<-data[unlist(map),]
  #Full Controls
  try(model<-brglm(sig_minor_war_mh~bin*amnest_init+truth_init+prisr_init
                   +powtran_init+as.factor(incompatibility)+poly(spline, 3, raw=T), data=boot.sample))
  if(is.na(model)) next
  if(table(is.na(model$coefficients))[1]<length(coef(model))) next
  if(model$converged==F) next
  effects<-coef(model)
  coef<-rbind(coef, effects)
  #cat("\r", i)
  i<-i+1
  model<-NA
}

#Extract coefficients and 95\% CI
coef<-na.omit(coef)
tab<-rbind(coef(base),apply(coef,2,sort)[c(25,976),])

#Replace with 90\% CI  for effect w/o substitutes
tab[2:3,3]<-sort(coef[,3],decreasing = F)[c(50,951)]

tab<-round(tab, 3)
tab<-t(tab)
colnames(tab)<-c("est", "lb", "ub")

#Order estimates
tab<-tab[c(3,2,nrow(tab),4:(nrow(tab)-1), 1),]

colnames(m1[[2]])[1:3]<-c("est", "lb",
                          "ub")
m1[[2]][1:3]<-round(m1[[2]][1:3], 3)
tab<-rbind(tab, m1[[2]][,1:3])

tab<-data.frame(tab)
tab$CI<-paste("[", tab$lb, " --- ", tab$ub, "]", sep="")

tab<-as.vector(t(tab[,c("est", "CI")]))

tab<-rbind(cbind(tab), cbind(round(rbind(nobs(base), length(unique(base$data$pam_caseid)), logLik(base)), 3)))
tab<-data.frame(tab)


rownames(tab)[c(seq(1,(nrow(tab)-3),2), (nrow(tab)-2):nrow(tab))]<-c("Amnesty Implementation", "Any Substitutes",
                 "Amnesty Implementation $\ times$ Any Substitutes",
                 "Truth Implementation",
                 "Prisoner Release Implementation",
                 "Power-sharing Implementation",
                 "Incompatibility", "Time Since Last War",
                 "Time$^2$", "Time$^3$", "Intercept",
                 "Marginal Effect of Substitutes without Amnesty",
                 "Marginal Effect of Substitutes with Amnesty",
                 "Difference", "N", "Clusters", "Log-Likelihood")
tab
names(tab)<-"Effect on War Recurrence"
tab<-print(xtable(tab, digits=3, align=c("l","c")), hline.after = c(25,28))
write(tab, file="Table 7.tex")

#With FORGE Substitutes - Appendix Figure 16 #####
forge_form<-as.formula(paste("sig_minor_war_mh~as.factor(bin_forge)*as.factor(amnest_init)+",
                             paste(controls, collapse="+")))
forge<-brglm(forge_form, data=data)
set.seed(1015)
forge_boot<-cluster_bootstrap(forge, n=1000, amn="amnest_init", bin="bin_forge")
m2<-int_plot(forge, forge_boot, amn="amnest_init", bin="bin_forge")
m2[[1]]

pdf("Fig16.pdf", width=7, height=4)
print(m2[[1]])
dev.off()

#With Substitutes Measured Annually - Appendix Figure 18 #####
set.seed(1015)
yearly_form<-as.formula(paste("sig_minor_war_mh~as.factor(bin_yearly)*as.factor(amnest_init)+",
                              paste(controls, collapse="+")))
yearly<-brglm(yearly_form, data=data)
yearly_boot<-cluster_bootstrap(yearly, n=1000, amn="amnest_init", bin="bin_yearly")
m3<-int_plot(yearly, yearly_boot, amn="amnest_init", bin="bin_yearly")

m3[[1]]
pdf("Fig18.pdf", width=7, height=4)
print(m3[[1]])
dev.off()

#With Substitutes without Infighting Condition - Appendix Figure 17#####
noinform<-as.formula(paste("sig_minor_war_mh~as.factor(bin_noinf)*as.factor(amnest_init)+",
                           paste(controls, collapse="+")))
noinf_mod<-brglm(noinform, data=data)
set.seed(1015)
noin<-cluster_bootstrap(noinf_mod, n=1000, bin="bin_noinf", amn="amnest_init")

m5<-int_plot(noinf_mod, boot.out=noin, amn="amnest_init", bin="bin_noinf")
m5[[1]]
pdf(file="Fig17.pdf", width=7, height=4)
print(m5[[1]])
dev.off()

#Including First Stage Vars - Appendix Figure 19#####
r3<-c("v2x_cspart", "v2x_polyarchy", "v2x_rule")
r3plus<-c(controls, r3)

R3_form<-as.formula(paste("sig_minor_war_mh~as.factor(bin)*as.factor(amnest_init)+",
                          paste(r3plus, collapse="+")))
R3_mod<-brglm(R3_form, data=data)
set.seed(1015)
R3_boot<-cluster_bootstrap(R3_mod, n=1000, bin="bin", amn="amnest_init")
R3_result<-int_plot(R3_mod, boot.out = R3_boot, amn="amnest_init", bin="bin")
R3_result[[1]]
pdf("Fig19.pdf", width=7, height=4)
print(R3_result[[1]])
dev.off()


#Including DDR Variables - Appendix Figure 20 #####
r1<-c("disarm_init","demob_init","reint_init")
r1plus<-c(controls, r1)
R1_form<-as.formula(paste("sig_minor_war_mh~as.factor(bin)*as.factor(amnest_init)+",
                          paste(r1plus, collapse="+")))
R1_mod<-brglm(R1_form, data=data)
set.seed(1015)
R1_boot<-cluster_bootstrap(R1_mod, n=1000, amn="amnest_init", bin="bin")
m4<-int_plot(R1_mod, forge_boot, amn="amnest_init", bin="bin")
m4[[1]]
pdf("Fig20.pdf", width=7, height=4)
print(m4[[1]])
dev.off()

#Bayesian Model Averaging - Appendix Figures 21-22 #####
bma_dat<-na.omit(data[,colnames(data) %in% c("sig_minor_war_mh", "bin", 
                                                       "amnest_init", r3plus, r1plus, "incompatibility")])
#colnames(bma_dat)[c(3,2)]<-c("Amnesty Implementation", "Any Substitutes")
colnames(bma_dat)[which(colnames(bma_dat)=="amnest_init")]<-"Amnesty Implementation"
colnames(bma_dat)[which(colnames(bma_dat)=="bin")]<-"Any Substitutes"

vars<-unique(c(r1plus, r3plus))
vars<-vars[!vars %in% c("spline", "spline2", "spline3")]
set.seed(1015)
bas<-bas.glm(as.formula(paste("sig_minor_war_mh~`Any Substitutes`*`Amnesty Implementation`+
                              as.factor(incompatibility)+poly(spline,3,raw=T)+",
                              paste(vars, collapse = "+"))),
             data=bma_dat, family="binomial"(link="logit"), betaprior = g.prior(3), force.heredity = T)
plot(coef(bas), subset=3, ask=F)

pdf("Fig21.pdf", width = 6, height = 3)
print(plot(coef(bas), subset=3, ask=F))
dev.off()

plot(coef(bas), subset=bas$n.vars, ask=F)

pdf("Fig22.pdf", width = 6, height = 3)
print(plot(coef(bas), subset=bas$n.vars, ask=F))
dev.off()

summary(bas)

rm(list=setdiff(ls(), c("cluster_bootstrap", "int_plot")))